Event-by-event hydrodynamics and elliptic flow from fluctuating initial state 



O 
(N 

C 
c3 



H. Holopainen a ^El H. Nicmi c and K.J. Eskola '*! 

a Department of Physics, P.O.Box 35, FIN-40014 University of Jyvaskyla, Finland 

b Helsinki Institute of Physics, P.O.Box 64, FIN-00014 University of Helsinki, Finland 

'Frankfurt Institute for Advanced Studies, Ruth-Moufang-Str. 1, D-60438 Frankfurt am Main, Germany 

(Dated: 3 January, 2011) 

We develop a framework for event-by-event ideal hydrodynamics to study the differential elliptic 
flow which is measured at different centralities in Au+Au collisions at Relativistic Heavy Ion Collider 
(RHIC). Fluctuating initial energy density profiles, which here are the event-by-event analogues of 
the eWN profiles, are created using a Monte Carlo Glauber model. Using the same event plane 
method for obtaining V2 as in the data analysis, we can reproduce both the measured centrality 
dependence and the pt shape of charged-particle elliptic flow up to pr ~ 2 GeV. We also consider 
the relation of elliptic flow to the initial state eccentricity using different reference planes, and discuss 
the correlation between the physical event plane and the initial participant plane. Our results 
demonstrate that event-by-event hydrodynamics with initial state fluctuations must be accounted 
for before a meaningful lower limit for viscosity can be obtained from elliptic flow data. 

PACS numbers: 25.75.-q, 25.75.Dw, 25.75.Ld, 47.75. +f 
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I. INTRODUCTION 

Azimuthal anisotropy of final state particles produced 
in ultrarelativistic heavy ion collisions can be used to 
measure the collective behavior of the dense particle 
system formed in such collisions [l(. The strong az- 
imuthal anisotropy, which has been observed in the trans- 
verse momentum spectra of hadrons in Au+Au collisions 
at the Relativistic Heavy Ion Collider (RHIC) of the 
Brookhaven National Laboratory, is also a signature of 
the formation of strongly interacting partonic matter, the 
Quark-Gluon Plasma (QGP). 

Ideal hydrodynamics has been successful in predicting 
and explaining the measured elliptic flow in Au+Au colli- 
sions at RHIC [2Hl3| . Currently, a lot of effort is devoted 
for developing a description of the QCD-matter evolu- 
tion in terms of dissipative hydrodynamics. The recent 
results show that even a small viscosity can considerably 
decrease the elliptic flow pdH19j . 

However, all these ideal and viscous hydrodynamic 
studies tend to underestimate the elliptic flow in most 
central collisions. Generally, the explanation for the 
deficit has been thought to be the initial state density 
fluctuations which have not been accounted for. In addi- 
tion to taking into account the density fluctuations them- 
selves, special care should be taken in computing the el- 
liptic flow with respect to the same reference plane as in 
the data analysis. 

The initial state fluctuations can be implemented e.g. 
via a Monte Carlo Glauber (MCG) model which makes 
possible to study the fluctuations of the initial matter 
eccentricity. Geometric fluctuations in the positions of 
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nucleons have been shown to increase the initial eccen- 
tricity, which is then suggested to translate into elliptic 
flow of final state particles |2fJ . Furthermore, the refer- 
ence plane plays a crucial role: the eccentricity is larger if 
one calculates it using the participant plane (determined 
by the transverse positions of the participant nucleons 
and the beam axis) instead of the reaction plane (deter- 
mined by the impact parameter and the beam axis). 

Recently, in Ref. [21(, hydrodynamical calculations 
were performed using averaged initial density profiles 
which were obtained from MCG calculations. Before av- 
eraging over the profiles, the transverse coordinate axes 
were rotated in each event so that the participant planes 
were on top of each other. In this manner it is possible to 
get an averaged initial profile that takes into account the 
eccentricity fluctuations in the initial state. For Au+Au 
collisions at RHIC, however, the effects of such plane ro- 
tations on the integrated V2 were small. 

While the above studies are steps to the right direction, 
it is obvious that without doing event- by- event hydrody- 
namic simulations, it is impossible to know how closely 
the computational participant plane corresponds to the 
physical event plane which is determined from the ob- 
served final state hadron momenta. 

So far, genuine event-by-event models where hydro- 
dynamics is run event by event using fluctuating initial 
density profiles, have been presented in Refs. [23 - l27| . In- 
terestingly, a similar two-particle correlation ridge as ob- 
served in the experiments [28[, is seen to form into the 
rapidity-azimuth- angle -plane both in NeXSphcrio [24j 
and more recently in Ref. [25(. This suggests that the 
puzzling ridge may well be another consequence of the 
fluctuations in the initial state. 

Also higher flow coefficients have been measured |29l - 
l3l| and recent studies [32| show that the initial state 
density fluctuations may play an important role in under- 
standing the centrality dependence of the ratio Vi/(v2) 2 - 
Triangular flow arising from event-by-cvent fluctuations 



|33| is also one of the things that should be studied fur- 
ther with event-by-event hydrodynamics. 

In this paper, we introduce an event-by-event ideal hy- 
drodynamics framework to study the following ^-related 
problems: With ideal hydrodynamics using averaged ini- 
tial states, (i) there is a vi deficit in central collisions, 
as discussed above; (ii) the shape and centrality depen- 
dence of i>2 (pt) are unsatisfactory in that the pr slopes of 
V2 easily become too steep and elliptic flow increases too 
much towards noncentral collisions; (in) elliptic flow is 
computed relative to the initial reaction plane or in the 
best case to the participant plane |21[ but not relative 
to the event plane, which is commonly used in the ex- 
periments; (iv) one does not know how closely the event 
plane and the initial participant plane correspond to each 
other. A concrete illustration of the problems (i)-(ii) can 
be found in Fig. 7.5. of Ref. I34|. and also in Fig. 5 of 
our previous elliptic flow study [12j . 

We will show how event- by- event ideal hydrodynamics, 
initiated with a fluctuating initial density profile obtained 
from a MC Glauber model, and especially the determi- 
nation of vi with respect to the event plane, conveniently 
solves the problem of the V2 deficit in the most central 
Au+Au collisions at RHIC. Simultaneously, we can sig- 
nificantly improve the agreement with the data for v-2 
at all centrality classes up to 30-40% most central colli- 
sions in the typical applicability region of hydrodynam- 
ics, pt < 2 GeV. This in turn has the very important im- 
plication that viscous effects can in fact be allowed to be 
smaller than previously thought. Finally, we also show 
the correspondence between the event and participant 
planes and study the relation between the elliptic flow 
and initial eccentricity using different reference planes. 

The rest of the paper is constructed as follows: First, 
in Sec. II we introduce our framework for event-by-event 
hydrodynamics. Details discussed there are our MCG 
model, computation of the fluctuating initial energy den- 
sity profiles, MC modeling of thermal spectra of final 
state hadrons, and MC modeling of the resonance de- 
cays. We also try to discuss the points where our model- 
ing could be improved. Section III is devoted for defining 
the event plane and elliptic flow. Also eccentricity issues 
are discussed there. Our results are presented in Sec. IV 
and conclusions are given in Sec. V. 



dius and d — 0.54 fm for the surface thickness. In the 
transverse (x, y) plane the two nuclei are separated by an 
impact parameter b between the centers of mass of the 
nuclei, which is determined by sampling the distribution 
dN/db oc b in the region < 6 < b max = 20 fm > 2ro. 
The longitudinal z coordinate is taken into account when 
sampling the initial nucleon positions but in what follows 
it does not play any role. 

Nucleons i and j from different nuclei are then assumed 
to collide if their transverse distance is small enough, 
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where <7nn is the inelastic nucleon-nucleon cross section. 
We apply here cjnn = 42 mb for Au+Au collisions at 
JsWn = 200 GeV. 

We note that our simple MCG model fails to reproduce 
the correlations between the nucleons, since we use the 
WS distribution for determining the nucleon positions 
independently from each other. In [35| it was observed 
that a realistic model, which accounts for nucleon corre- 
lations [36| , can be well approximated using an exclusion 
radius which prevents nucleon overlap. Using such ra- 
dius, or giving a finite size for the nucleons [2l| . causes 
deviations from the WS distribution which should then 
be compensated by tuning of the parameters in the ini- 
tially sampled WS distribution. 

To keep our modeling as transparent as possible we, 
however, choose not to apply an exclusion radius or a nu- 
cleon size in our MCG model since according to Ref. [21| 
only a f 0% uncertainty in the initial eccentricity can be 
expected, which is a much smaller effect than e.g. the 
overall uncertainties related to the choice of the initial 
density profiles. 

Next, we define the centrality classes using the num- 
ber of participant nucleons, A part , for simplicity. We 
have plotted the distribution of events as a function of 
-^Vpart in Fig. [TJ As indicated there, we slice our total 
event distribution in N-oart so that each iV Dart interval 



*part 



part 



corresponds to a centrality class which contains a certain 
percent of total events. The impact parameter may thus 
freely fluctuate within each centrality class. 



II. EVENT-BY-EVENT HYDRODYNAMICS 
FRAMEWORK 



A. MC Glauber model and centrality classes 

We use here a MCG model to define the centrality 
classes and to form initial states with fluctuating density 
profiles. First, we distribute the nucleons in the colliding 
nuclei randomly using the standard, spherically symmet- 
ric, two-parameter Woods-Saxon (WS) nuclear density 
profile as the probability distribution. Our WS param- 
eters for the gold nucleus are Ra ~ 6.37 fm for the ra- 



B. Initial density profiles 

In order to utilize the MCG-given initial state to start 
hydrodynamics, we must next somehow transform the 
positions of the wounded nucleons or binary collisions 
into energy density or entropy density. These would 
be the fluctuating event- by- event MCG analogues of the 
conventional eWN, eBC and sWN, sBC Q average ini- 
tial densities. For simplicity, we consider here just the 
cWN-type of profile and leave the profile fine tuning for 
future work. The energy density is now distributed in 
the (x, y) plane around the wounded nucleons using a 2D 
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FIG. 1: Our definition of centrality classes for Au+Au col- 
lisions at y / ijviv = 200 GeV. Distribution of the number of 
participants is calculated from a Monte Carlo Glauber model 
without a nucleon exclusion radius. 



Gaussian as a smearing function, 
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where K is a fixed overall normalization constant and a 
is a free smearing parameter controlling the width of our 
Gaussian. In each event, the impact parameter defines 
the direction of the x axis and the origin of the (x, y) 
plane is determined so that the energy-density weighted 
coordinate averages become (x) = (y) = fm. 

For the hydrodynamical description to be meaningful, 
the initial state should not have too sharp density peaks. 
In our MCG model we have given an effective interaction 
radius yAAAv/Tr/2 ~ 0.6 fm for the colliding nucleons, 
which sets a natural order of magnitude for a. To probe 
the sensitivity of our results to the initial state smear- 
ing, we will consider two values, a = 0.4 fm and 0.8 fm. 
With the current setup, we cannot reduce a further, as 
this would require a smaller step size in our hydrodynam- 
ical code, and consequently much more CPU time. One 
should then also develop a way to handle multiple sep- 
arate freeze-out surfaces, see the discussion in Sec. Ill CI 
These developments we leave as future improvements. 

The reason to choose the energy density to be smeared 
rather than the entropy density, is mostly technical and 
due to the fact that our focus here is on understanding 
the transverse flow phenomena. Since we now avoid us- 
ing the Equation of state in forming the initial energy 



density profiles in each event, we have a more direct con- 
trol on the input energy density (pressure) gradients that 
drive the evolution of the transverse flow and its asym- 
metries. In our case, the total energy per rapidity unit in 
each event, J dxdy e(x,y), thus remains independent of 
a, while the total entropy per rapidity unit and thereby 
also the final state multiplicity depend on a . 

For Au+Au collisions at ^/smn — 200 GcV, we use 
the value K = 37.8 GeV/ fm. With this, we reproduce 
the initial total entropy of Ref. [12| when averaging over 
many initial states in central (b = 0) collisions when 
a = 0.4 fm. Motivated by the EKRT minijet (final state) 
saturation model 37J and Ref. [12| , we fix the initial time 
to ro = 0.17 fm for all events. 



C. Hydrodynamics, freeze out and resonance 
decays 

For obtaining the ideal-fluid hydrodynamic evolution 
of the system, we solve the standard equations 



d»T» v 
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together with an Equation of State (EoS) which re- 
lates pressure with the energy density and net-baryon 
number density, P = P(e,ns)- As we are interested 
in particle production at mid-rapidity, we assume the 
net-baryon density to be negligible. Since the rapid- 
ity distributions of hadrons are approximately flat at 
mid-rapidities we can safely simplify our hydrodynami- 
cal equations by assuming longitudinal boost-invariancc. 
We solve this (2+l)-dimensional numerical problem us- 
ing the SHASTA algorithm [H [H[ which is also able to 
handle shock waves. 

As the Equation of State (EoS), we choose the EoS 
from Laine and Schroder [4fJ] . At high temperatures this 
EoS has been matched with the lattice-QCD data and at 
low temperatures with a hadron resonance gas containing 
particles of mass m < 2 GcV. This EoS has a "cross- 
over" transition from the QGP to the hadron gas. 

Thermal spectra for hadrons are calculated using the 
conventional Cooper- Frye method [4JJ, where particle 
emission from a constant-temperature surface a is cal- 
culated according to 



dN 
d 2 prdy 



= / f{x,p)p' 1 da fJ 



(4) 



where f(x,p) is the particle number-distribution func- 
tion in momentum at a certain space-time location. The 
frcczc-out temperature Td oc = 160 MeV is fixed so that 
we reproduce the measured pt spectrum of pions [42j 
when averaged initial states are considered. 

Our surface finding algorithm operates in the (r, r)- 
plane for all spatial azimuthal angles. Currently, we can 
find only surfaces which go through r = 0. Due to the 
initial state fluctuations there might simultaneously exist 
also other, disconnected, freeze-out surfaces which our al- 
gorithm does not recognize. Wc have checked that for the 



centrality classes and smearing parameters a considered 
here, only a a few percent of the events actually con- 
tain such a surface. In any case, since these additional 
surfaces typically originate from a few-nucleon collisions, 
they contribute negligibly to particle production in not 
too peripheral Au+Au collisions. Making a smaller can 
also increase the number of disconnected freeze-out sur- 
faces. To ensure the applicability of our framework, we 
prefer not to consider centrality classes more peripheral 
than 30-40% or a < 0.4 fm in the present study. 

For the flow analysis, we need individual final state 
particles. In generating these using the computed ther- 
mal spectrum as the probability distribution, we assume 
the total number of thermal particles in a rapidity unit to 
be fixed individually in each event. The transverse mo- 
mentum (p x ,p y ) for each particle is thus sampled from 
the distribution dN/cPprdy calculated in Eq. (QJ. Due 
to the assumed boost-symmetry, we are not equipped to 
consider rapidity distributions, thus y is sampled from a 
flat distribution in the interval \y\ < 0.5. 

Note that above we have neglected the fluctuations in 
the number of emitted thermal particles. In principle 
one could derive these fluctuations separately from the 
thermal distributions for each freeze-out surface element. 
However, it is not so clear how to treat the space-like 
parts of the surface in this case. Since in the collisions 
considered here there are of the order of 1000 particles 
per unit rapidity, these fluctuations can in any case be 
expected to be negligible in comparison with the initial 
state fluctuations. 

Once we have generated all the thermal hadrons, we 
still need take into account the strong and electromag- 
netic decays. We let the thermal resonances decay one 
by one using PYTHIA 6.4 [43[. Some decay products 
can fall outside our rapidity interval |y| < 0.5. On the 
other hand, there would also be decay products arriving 
from \y\ > 0.5 which we do not consider here. We have 
checked that instead of increasing the width of our ther- 
mal particle rapidity window, to speed up the analysis, 
we can simply count all decay products into our rapidity 
acceptance regardless of their actual rapidity. 



D. Event statistics 

Our main goal is to compare the event- by- event hy- 
drodynamic results with the ones obtained by more con- 
ventional non-fluctuating hydrodynamics initiated with 
averaged initial states. 

For event-by event hydrodynamics, we make 500 hydro 
runs in each centrality class. This amount of hydro runs 
seems enough for the hadron spectra and elliptic flow 
analysis. To increase statistics we make 20 final state 
events from every hydro run, thus we have 10 000 events 
in total. To check that using each hydro run 20 times is 
sensible, we have checked that doing 250 hydro runs and 
40 events from each leads to the same flow results. 

To create an averaged initial state, we sum together 



20 000 initial states generated by our MCG model. Such 
large number of events is required since fluctuations near 
the edges of the system easily affect the final value of 
elliptic flow if the density profile is otherwise smooth. We 
then do one hydro run with the averaged initial state for 
each centrality class. To make a fare comparison with the 
event-by-event hydro results, we do the resonance decays 
and analysis using the same code for the averaged initial 
state case as for the event-by-event hydro case, making 
10 000 final state events from this one hydro run. 



III. ELLIPTIC FLOW ANALYSIS 

A. Elliptic flow and event plane 

The transverse momentum spectra of hadrons can be 
written as a Fourier series, 
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where <f> is the hadron momentum's azimuthal angle with 
respect to the reaction plane defined by the impact pa- 
rameter. The flow coefficients v n can then be computed 
from 
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When we have fluctuations in the initial state, calcula- 
tion of v n is not so straightforward. In the hydrodynamic 
runs, where we always know the direction of our impact 
parameter, we can calculate the elliptic flow with respect 
to the reaction plane. If we want to compare with ex- 
periments, we should use the same analysis methods and 
definitions as in the data analysis. In this work we use 
the event plane method [44|, |45[ which is a common way 
to calculate V2- Since it is not (yet) typically used in hy- 
drodynamical calculations, let us briefly recapitulate the 
main points (see Ref. |44| for details). 

We first define an event flow vector Q n for the nth 
harmonic. The event flow vector in the transverse plane 
is 



Q, 






{pTi cos(n<j>i),pTi sm(n<t>i)), 



(7) 



where we sum over every particle in the event and where 
4> is measured from the x axis which is here fixed by the 
impact parameter. The event plane angle ip n for each 
event is then defined to be 



''/•„ 



<vccitXn[CJ n y/CJ niX ) 



(8) 



with arctan placed into the correct quadrant. The "ob- 
served" v n is calculated with respect to the event planes 
obtained above, 



u„{obs} = ((cos(n(oX - ^n))))€ 



(9) 



where the inner angle brackets denote an average over 
all particles i in one event and the outer ones an average 
over all events. In order to remove autocorrelations, the 
particle i is excluded from the determination of the event 
flow vector when correlating it with the event plane. 

Since in our finite rapidity interval we have only a finite 
number of particles available for the event plane deter- 
mination, the obtained event plane fluctuates from the 
" true" event plane. (In our event- by- event hydrodynam- 
ics, the true event plane in each event would correspond 
to the average event plane obtained by generating in- 
finitely many final states from one hydro run.) The ob- 
tained v„{obs} is corrected using the event plane resolu- 
tion for the harmonic n 



n„ 



(cos(n(^-^ rue ))>. 



(10) 



where tp^ UG defines the true event plane and the an- 
gle brackets stand for an average over a large sample of 
events. Because experimentally it is not possible to find 
the true event plane, the event plane resolution must be 
estimated. 

In the two-subevents method, which also we will use, 
each event is randomly divided into two equal subevents 
A and B. The event plane resolution for each of these 
subevents is then \4M 



n: 



yj (casing -il>B))). 



(11) 



If the fluctuations from the true event plane are Gaussian, 
one can analytically obtain the following result [44j 



B. Initial eccentricity and participant plane 

The reaction plane eccentricity of the hydrodynamical 
initial state can be defined as (see e.g. Ref. [201 ] ) 



(14) 



where 
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where the averaging is done over the energy density pro- 
file of Eq. ©. 

Since the positions of wounded nucleons, however, fluc- 
tuate from one event to another, tilting the transverse 
coordinate axes suitably we can actually get a larger ec- 
centricity than epp above. Thus it is not so clear what 
the most correct reference plane should be. 

The reference plane that maximizes the initial eccen- 
tricity can be expected to correlate better with the event 
plane than the reaction plane. For this purpose, one may 
define the participant eccentricity [20J 
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where a xy = (xy) — (x) (y) . In this case the reference 
plane is the participant plane which is defined by the z 
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X„exp(- X 2/ 4 ) J ( x 2 /4)+Jl(x 2 /4) ( (12) 



where Iq and I\ are modified Bessel functions and Xn ~ 
y/~N, with N referring to the number of particles. Since 
we can calculate TZ^ h from the subevents, we can nu- 
merically solve Xn" b horn Eq. (fi"2|) . Because the number 
of particles in the subevents is half of those in the full 
events, x^ 1 " — V%X™ h > an d we can calculate the resolu- 
tion 7?4 u11 for the full events. Finally, the flow coefficients 
are obtained as 



f„{obs} 
■R£ li ' 



(13) 



The elliptic flow results computed with this method 
are denoted here as i>2{EP}. We also compute the ellip- 
tic flow from Eq. © with respect to the reaction plane 
using both fluctuating and averaged initial states. In 
the reaction plane case we have no corrections coming 
from statistical fluctuations. These results are denoted 
as 1>2{RP} in what follows. 
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FIG. 2: (Color online) Transverse momentum spectra of pos- 
itive pions for Au+Au collisions at ^/snn = 200 GeV cal- 
culated with averaged and fluctuating initial states varying 
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axis (beam direction) and the x axis which is first rotated 
around the z axis by the angle 
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In what follows, we will compute the elliptic flow also 
with respect to the participant plane, 



W 2 {PP} = ((C0S(2(& - Vpp))))cvonts 



(18) 



and consider the relation of elliptic flow to the initial 
eccentricity using both the reaction plane and the par- 
ticipant plane as the reference. 



IV. RESULTS 

Below, we present the results for pion spectra, elliptic 
flow, eccentricities and the correlation of the event and 
participant planes. The genuine event- by- event calcula- 
tions using smearings a = 0.4 and 0.8 fm, are compared 
with the results obtained using an averaged initial state. 

In Fig. [2] we show the pt spectra of positive pions 
from these three different hydro calculations and from 
the PHENIX collaboration [42[. As explained in Sec. Ill Bl 
our multiplicity depends on the Gaussian smearing width 
a, hence the (small) difference between the points with 
a = 0.4 fm and a = 0.8 fm at low px- 

We can also see that at higher pt we get more particles 
with the fluctuating initial states than with the averaged 
initial state case. This follows from the fact that in the 
fluctuating initial states there are larger pressure gradi- 



ents present. For the same reason, the high-px spectra 
are quite sensitive to the value of a: with a larger a, the 
pressure gradients are smaller and the pr spectra steeper. 
This is in fact an interesting observation, suggesting that 
with fluctuating initial states the applicability region of 
(event- by- event) hydrodynamics may extend to higher pt 
than previously (see e.g. Refs. 4f||47|) thought. In any 
case, the obtained pt spectra agree with the data suffi- 
ciently well, so that we can meaningfully next study the 
elliptic flow. 

In Fig. [3] we plot the elliptic flow of charged particles 
as a function of pr at different centralities. We show the 
event- by- event results for U2{EP}, U2{RP} and U2{PP}, 
as well as i>2{RP,ave} which is obtained from averaged 
initial states. 

First, we observe, that i>2{RP} and i?2{RP,ave} are 
quite close to each other (although in the panel (c) some 
statistical fluctuations seem to be still present), and cspc- 
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FIG. 4: (Color online) Elliptic flow of charged particles as 
a function of pr at different centralities with two different 
values for Gaussian smearing parameter a. 



cially that in central collisions there is a significant deficit 
of V2 relative to the data. Second, we see that our U2{EP} 
agrees very well (within the estimated errors) with the 
data up to pr ~ 2 GeV in all centrality classes. No- 
tice also the difference between the uncorrected «2{obs} 
and the corrected, final, V2{EP}; especially for central 
collisions, the corrections 7?. 2 are quite large. Thus, fluc- 
tuations alone are not sufficient in explaining the vi data 
but that in addition to taking into account the fluctu- 
ations - the computed V2 must be defined in the same 
way as in the experimental analysis. 

Third, we notice that the relative increase from 
i'2{RP} to «2{EP} decreases from central to peripheral 
collisions: w 2 {EP}/w 2 {RP} = O(10) in panel (a) and 
0(1.2) in panel (d). Fourth, contrary to our original ex- 
pectation, vi {RP,ave} for semi-peripheral collisions is ac- 
tually below (and not above) the data at pr ~ 1.5 GeV. 
This is due to the fact that with our MCG model and 
smearing, the actual energy density profiles become flat- 
ter and less eccentric than the conventional eWN profiles 
obtained from an optical Glauber model. As a result, we 
get a smaller i^lRP} than e.g. in Ref. [13], and thus also 
in the 20-30% centrality class there is room for an in- 
crease from «2{RP} to i>2{EP}. From these observations 
we can conclude that we have answered the questions 
(i)-(iii) presented in Sec. [U 

Fourth, Fig. [^indicates that W2{PP} is very close to 
w 2 {EP} in all centrality classes. This result suggests that 
the participant plane indeed is quite a good approxima- 
tion for the event plane. 

In Fig. |4] we show the effects of varying our Gaussian 
smearing parameter a. We see that our elliptic flow re- 
sults are quite insensitive to a: Doubling the value of a 
causes only of the order of 10% changes in our v^^Pt)- 
We remind, however, that our p^-spectra and multiplic- 
ity of pions were not as stable against a but we expect 
that doing more proper fitting to the pion spectra by fine- 
tuning Tdec and the initial overall normalization constant 
K, would not affect our vi results significantly. 

In Fig. [5]we plot the integrated elliptic flow for the four 
different cases considered above and the data from the 
PHOBOS collaboration [50(. As expected on the basis 
of Figs. [3] and [2 our results i>2{EP} and i^jPP} now 
agree with the data very well, while the i^jRP} results 
fall significantly below the data. 

Next, Fig. [5] shows the computational quantity v-2.lt 
which is often discussed. In the PHOBOS result [5fJ], V2 
is determined relative to the event plane while the initial 
state eccentricity is computed relative to the participant 
plane. We reproduce the PHOBOS W2/e if we do the 
same, i.e. use epp from Eq. (fl~6|) . Interestingly, if wc 
replace both the elliptic flow and the eccentricity by their 
reaction plane analogues, we can still get a scaling law 
that agrees with our w 2 {EP}/epp and with the data. This 
figure illustrates again the importance of the consistency 
in the reference plane definition. 

Finally, we answer the question (iv) presented in Sec. [I] 
Figure [7] shows the correlation between the event plane 



and the participant plane as well as the correlation be- 
tween the event plane and the reaction plane. We plot 
the distribution of events as a function of the angle dif- 
ferences -02 — ippp and -02 — ipRP- For this figure, we 
have used each hydro run only once. We notice that 
in central collisions, the planes are more weakly cor- 
related than in semi-peripheral collisions where clearer 
peaks around tp2 = ippp, V'rp arise. As expected based on 
Fig- 13 the participant plane is indeed found to be quite a 
good approximation for the event plane, in all centrality 
classes. However, fluctuations of the event plane around 
the "true" event plane are much larger in central colli- 
sions and thus the correlation between the event plane 
and the participant plane in Fig. [7] looks weaker for cen- 
tral collisions than for the more peripheral ones. 



V. CONCLUSIONS 

The main result of this paper is that using event-by- 
event ideal hydrodynamics with MCG-generated fluctu- 
ating initial density profiles we can simultaneously re- 
produce the measured centrality dependence and the px 
shape of charged-particlc elliptic flow up to pt ~ 2 GeV. 
Also the measured pion spectra are quite well repro- 
duced, although we have not made an effort to fine- 
tune the model parameters. In particular, in addition 
to accounting for the fluctuations in the system, we have 
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FIG. 6: (Color online) Integrated elliptic flow of charged par- 
ticles for Au+Au collisions at ^/sjvjv = 200 GeV divided by 
initial eccentricity. Theoretical calculations correspond to 
fluctuating initial states with a = 0.4 fm. Data from the 
PHOBOS collaboration [50] are shown with statistical and 
systematic errors added in quadrature. 
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averaged initial states are shown for a — 0.4 fm. Data from 
the PHOBOS collaboration [50] are shown with statistical and 
systematic errors added in quadrature. 



demonstrated the importance of using the same V2 defi- 
nition as in the data analysis. 

We have performed all hydrodynamic simulations with 
zero viscosity. Thus, our results suggest that extracting a 
non-zero lower limit for the viscous coefficients from the 
measured V2{pr) of charged hadrons is practically impos- 
sible without further constraints to the model, especially 
to the initial state. We would like to emphasize, that 
we have for simplicity considered only the event- by- event 
analogues of the eWN initial profiles whose eccentrici- 
ties are typically smaller than e.g. those of the eBC- or 
CGC-type |5ll - [53| of profiles. Whether the data are still 
consistent with non-zero viscosity with these initial con- 
ditions, is left as a future exploration. Nevertheless, our 
results demonstrate that event- by- event hydrodynamics 
with initial state fluctuations must be accounted for be- 
fore a more reliable lower limit for viscosity can be ob- 
tained from elliptic flow data. 

We have shown that the definition of the reference 
plane with respect to which one determines V2, plays an 
important role especially in central collisions. On the one 
hand, if V2 is computed relative to the reaction plane (de- 
termined by the impact parameter), the fluctuating and 
averaged initial states lead practically to the same results. 
In this sense, the previous conventional ideal hydrody- 
namical results for the system evolution are still relevant 
in central enough collisions but one should not compare 
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decoupling temperature may vary from event to event. 
Instead of a fixed T^ ec applied here, one could imple- 
ment a dynamical freeze-out criterion as was done e.g. 
in Ref. [541] . However, in order to improve upon the well- 
known problem of the proton px spectra when partial 
chemical equilibrium is not applied, one could couple our 
hydro to a hadron cascade afterburner which would han- 
dle also the resonance decays [2a - l27j |. Related to the 
initial state, one should more closely inspect the uncer- 
tainties due the assumed energy density smearing, which 
is an avoidable issue with event- by- event hydrodynam- 
ics. Here we found out that v 2 remained fairly insensitive 
to Gaussian smearing width while pion px spectra were 
more sensitive to it towards larger pr- Also other pos- 
sible smearing functions should be studied. One should 
also consider a dynamical QCD-based model for the ini- 
tial fluctuations, in which case also the absolute initial 
density profiles should be computable. These tasks, and 
considering the effects of fluctuations on other observ- 
ables, we leave as future developments. 
Acknowledgments 



FIG. 7: (Color online) Correlation of the event plane with 
the participant plane, and with the reaction plane at different 
ccntralities and with a — 0.4 fm. The lines are to guide the 
eye. 



the reaction-plane v 2 to the event-plane v 2 quoted by the 
experiments. On the other hand, according to our re- 
sults, the initial participant plane seems to be quite a 
good approximation for the event plane in the presence 
of hydrodynamically evolving density fluctuations. 

The present work can obviously be improved in many 
ways. Especially in event- by- event hydrodynamics the 
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